arXiv:1508.05330v2 [math.NA] 19 Nov 2015 


Variants of the Empirical Interpolation Method: symmetric 
formnlation, choice of norms and rectangnlar extension 

Fabien Casenave^’^, Alexandre Ern^ and Tony Lelievre^ 

^ IGN LAREG, Univ Paris Diderot, Sorbonne Paris Gite, 

5 rne Thomas Mann, 75205 Paris Gedex 13, France 
^ cnrrently at SafranTech, Rne des Jennes Bois, 

Ghateanfort, GS 80112, 78772 Magny-Les-Hameanx, France 
3 Universite Paris-Est, GERMIGS (ENPG), 

77455 Marne la Vallee Gedex 2, France 

Abstract 

The Empirical Interpolation Method (EIM) is a greedy procedure that constructs approxi¬ 
mate representations of two-variable functions in separated form. In its classical presentation, 
the two variables play a non-symmetric role. In this work, we give an equivalent definition 
of the EIM approximation, in which the two variables play symmetric roles. Then, we give 
a proof for the existence of this approximation, and extend it up to the convergence of the 
EIM, and for any norm chosen to compute the error in the greedy step. Einally, we introduce 
a way to compute a separated representation in the case where the number of selected values 
is different for each variable. In the case of a physical field measured by sensors, this is useful 
to discard a broken sensor while keeping the information provided by the associated selected 
field. 
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1 Introduction 

Consider a function f : X xy ^ M.. The Empirical Interpolation Method (EIM) [Ull] is an 
offline/online procedure that provides an approximate representation Id{f) of / in separated 
form, where the integer d denotes the number of terms in the representation. The offline stage 
of the EIM consists in selecting some points {xi,... ,Xd) G and {yi,... ,yd) G V in a 
greedy fashion, such that 


Xk+i = argmaxll (/ - 4(/)) (x, •)||Lo°(y), 

xG?i 

yk+i = argmaxi (/ - 4(/)) {xk+i,y)\, 


yay 


(la) 

(lb) 
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where Ikif) denotes the separated representation constructed with the k first points {xi,yi), 

I G {1, k}. The method is efficient when the error / — Idif) is small for reasonably small 
values of d. Some functions qi{y), I € ,(i} and a matrix B of size dx d depending on 

these points are also constructed, see [1] and [1] for details. The separated representation is 
then obtained as 

ld{f)ix,y)= Mx)qiiy), (2) 

l</<d 

where the functions Xi{x), I € {1, ...,d}, solve the linear system Ylm=i Bi,mXm{x) = f{x,yi), 

I € ,(i}. The function Idif) satishes the following interpolation property [H Lemma 1]: 

for all m € {1, d}, 

Idif)ix,ym) = fix,ym), for all X G T”, (3a) 

Idif)ixm,y) = fixm, y), for all y €y. (3b) 

In practice, the size d is not chosen a priori, and the greedy procedure stops when | (/ — Ikif)) ixk+i,yk+i)\ 
is small enough. Define U := {fix, ■),x G ff}. Elements of U are functions from y to M. 

Theorem 1.1 (Existence of the decomposition, [4|, Theorem 1]) Assume that the in¬ 
terpolation points are chosen according to P^- pT]) and that d < dimspan(ZY). Then, the 
separated representation m is well-defined. 

In 12], it is observed that Idif) oan be rewritten in an algebraically equivalent form as 

Idif)ix,y)= Di,mfixi,y)fix,ym), ( 4 ) 

where the matrix D depends on the points xi, ym, l,m ^ {1, ...,d}, and can be constructed 
recursively during the offline stage of the EIM. It is easy to check that (I3a1) - ()3bp is satished 
if and only \i D = F~'^, where = fixi,ym), which motivates an alternative presentation 
of the EIM based on Equation Q where the variables x and y play symmetric roles. The 
double summation in (jH) emphasizes this symmetric role; note that, for instance, Idif)ix, y) = 
Y.i<i<dMx)miy), with Xiix) = Y.i<m<dDi,mfix,ym) and qiiy) = fixi,y). 

Let ll-lly be a norm on y and suppose that the interpolation points are now selected as 

Xk+i = argmax||(/-4(/))(x,-)b, (5a) 

xGA, 

yk+i = argmaxi (/ - 4(/)) ixk+i,y)\, (5b) 

y&y 

the difference with (|lap - pbl) being the arbitrary choice for the norm H-Hy, in the first line, 
instead of just IHlLt’oCV)- Oiie can exchange the roles of x and y in the previous algorithm, 
leading to 


yk+i = argmax II (/-4(/)) (•,?/)IIa-, (6a) 

y^y 

Xk+i = argmax|(/-4(/))(x,yfc+i)|, (6b) 

xGA. 

for an arbitrary norm H-H;!' on X. In general, the couple ixk+i,yk+i) resulting from ([6aP - ()6bl) 
differs from the one obtained with (|5ap - (HE]). Choosing the L°°-norm in y and X in (|5al) 
and (l6aP respectively, we infer that 

(xfc+i,?/fc+i) = arg max | (/- 4(/)) (x, ?/)|, 

(x,y)(^Xxy 


2 


thus recovering the choice made in For this specific choice of L°°-norms on df and 

3^, (f5a|) - (|5bh is actually equivalent to (f6all - (|6bp . 

The hrst contribution of this work is the following result: 

Theorem 1.2 Assume that the interpolation points are chosen according to (|5ap - ()5bp or (I6ap - 
(160) and that d < dimspan(Z/f). Then, the separated representation (jH) with D = F ^, where 
Fi,m = f{xhVm), is well-defined and satisfies (l3aP - (l3bp . 

Theorem [Q extends Theorem o in two ways. First, it allows for a more general selection 
of the interpolation points. Second, the existence of the decomposition is ensured up to 
convergence. Since dimspan(W) is usually unknown, a practical consequence is that for all 
k € N, either f{x,y) = Ik{f){x,y) for all {x,y) € T” x T, or the greedy procedure can be 
pursued. Moreover, we observe that the proof of Theorem 11.21 given below is straightforward 
and does not rely explicitly on the Kolmogorov n-witdh. 

The second contribution of this work is that starting from (jl]), we devise a rectangular ex¬ 
tension of EIM. An interesting application of this extension is an improved recovery procedure 
if it is decided a posteriori that the values of / at some points yi are not reliable and should 
be discarded. This situation is motivated for instance by sensor failure in the context of the 
Generalized Empirical Interpolation Method (GEIM) [3], which we address at the end of this 
work. Note that discarding a posteriori some interpolation point yi is not straightforward in 
the standard EIM setting since the whole construction relies on couples of points and functions 
defined by induction. In the present setting, the interpolation points xi can be kept even if 
the points yi are discarded, thereby leading to a rectangular formulation. The key idea is to 
replace the matrix inversion by a pseudo-inverse to compute the rectangular matrix D to be 
used in (jH). 

2 Symmetric formulation of EIM and proof of Theorem 11.2 

2.1 Symmetric formulation 

Consider some interpolation points (xi,..., xfi) G and (yi,..., yfi) € and dehne the 
matrix F G such that Fi^rn = f{xi,ym), for all l,m £ {1,..., d}. 

Lemma 2.1 (Existence) If the matrix F is invertible, there exists a function Ifif) ■ X x 
y —> M m the form ([3|) satisfying the interpolation property (f3al) - (l3a1) . 


Proof We observe that for all mo G {1, and all x £ X, 


h{f){x,ymo) 


'y ^ -Dl,rnf {xi^yrno)f {x,yrn) — ^ ^ 

l<l,m<d l<m<d 


y ^ kll,mf {xi,ymo) 
l<l<d 


fix,ym)- 


Hence, if F^D = Id, (l3al) holds. Similarly, if DF^ = Id, (I3bp holds. Since F is invertible, we 
see that we can choose D = F~'^. □ 

Note that the EIM interpolation is thus given by Id{f){x, y) = J2i<i,m<diP~^)m,ifixi,y)f{x, ym)- 


2.2 Proof of Theorem 11.21 

We begin with the following result on the stopping criterion of the greedy procedure. 
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Lemma 2.2 (Stopping criterion) Consider a set of points {xi,ym)i<i,m<k+i o,nd assume 
that the matrix := {f{xi,ym))i<i,m<k is invertible. If 

if - h{f)) (xfc+i, ?/fc+i) / 0, 

then the matrix := {f{xi,ym))i<i,m<k+i is also invertible. 


Proof Assume that F^^^ is not invertible. This means that there exist (Ai,..., A^+i) not 
all zero such that, V/ € {1, , m f{xi,ym) = 0 , which writes 

^ >^raf{xi,ym) =-Xk+if{xi,yk+i). (7) 

l<7n<k 

Note that necessarily, A^+i 7 ^ 0, since otherwise the matrix F^ would not be invertible. Let 
us introduce the matrix = [F^y^ so that Ikif){x,y) = Y.i<i,m<k 
Let mo G {1,..., k}. From ([7]), we infer that 

^tmo Y1 ^rnfixi,ym) = ->^k+l Y1 ?/fc+l)' 

and since, for all m,mo € {!,...,/c}, exchanging summations 

in the left-hand side leads to 

^mo = “Afc_|_i Dij^^f(xi,yk+i). 

l<l<k 

We deduce that 

X] ^mof{x,ymo) = -h+l'^ Dtmofy^yk+l)fix,ymo) = -h+lIk{f){x,yk+l)- 

l<mo<k l<l<k l<7no<k 

Taking x = Xk+i and using again ([7]), we infer that 

'^A:+l/(^A:+l) l/fc+l) ^k+lkkif'ji.Xk+ljUk+l) 

and recalling that A^+i 7 ^ 0 , we conclude that 

f {Xk+ljUk+l) — 7fc(/)(xfc-|_i, ?/fc-|_i). 


□ 

It is clear from (f5a)l - (|5bp or (f6a)l - (|6bp that for all k < dimspan(7/), (/ — Ikif)) {xk+i,yk+i) 7 ^ 
0. Then, from Lemma 12.21 Ifc+i is well-defined, and Theorem O is proved. 

Corollary 2.3 (Termination) If {xk+i,yk+i) is selected according to (ISa)l - (f5bl) or (IUa| - (f6bl) . 
then I (/ - hif)) {xk+i,yk+i)\ = 0, implies that (/ - 4(/)) {x,y) = 0 for all {x,y) € T x T- 


Proof Suppose that (xfc+i, j/fc+i) is selected according to (l5a|)-([5b]). Then I (/ —/fc(/)) (xfc+i, ?/fc+i)| = 0 
implies by definition of yk+i that for all y & y, | (/ — Ikif)) ixk+i,y)\ = 0. This in turn 
means that \\ if — Ikif)) (xk+i,')\\y — 0- Then, by definition of Xk+i, for all x G T, 

II (/ — /*;(/)) (x, •)|| 3 ; = 0 , and thus / = Ikif)- The proof is similar in the case where 
ixk+i,yk+i) is selected according to (f 6 a|) - (l 6 bp . □ 

The practical consequence of the above results is that, for all /c G N, either | (/ — Ikif)) ixk+i,yk+i)\ = 0 
and / = Ikif) iCorollary 12.311 . or Ik+i can be constructed (Lemma 12 . 21 ) and the greedy pro¬ 
cedure can be pursued. 
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Remark 1 (Goal-oriented approximation) In a reduced-basis context, the EIM procedure is 
typically used to obtain a separated representation of the coefficients of a parametrized par¬ 
tial differential equation. As an example, let us consider the right-hand side f{x,y) of a 
PDE discretized using a Galerkin basis {(pi,... . Here, x is the parameter and y the 

space variable. The discretized right-hand side is then b{x)j = f{x,y)pj{y) dy and us¬ 
ing the EIM approximation, it actually can be written as a linear combination of functions 
of x: {Id{b){x))j = J^Id{f){x,y)pj{y)dy = ^m)/ q/( a^z, c??/- This 

separated representation is very important for the overall efficiency of the greedy procedure 
(see W)- Gonsider a finite-dimensional subspace of a Hilbert space with basis {pi, . 

Define b{x)j = f{x,y)ipj{y)dy for all j € and its EIM-based approximation 

{Id{b){x))j = J^Id{f){x,y)(pj{y)dy. If the EIM approximation Id has been constructed us¬ 
ing (f5a1) - (l5bp . the error in the greedy procedure is assessed on the quality of the approximation 
of f{x,y). A goal-oriented alternative is to assess the quality of the approximation of f on the 
object b to be approximated: 


Xk+i = argmax || (6 - Id{b)) (x)||, 
x^X 

yk+i = argmax] (/ - 4(/)) {xk+uy)\, 
y&y 


(8a) 

(8b) 


where jl-jl can be any norm on . 

3 Extension to rectangular cases 

3.1 Different number of interpolation points x and y 

Consider a given EIM approximation, where d couples {xm,ym), m £ Af := {1, ...,d} have 
been selected. Suppose now that for some reason, some points ymo^ "Iq € A/q C A/", must 
be discarded. For instance, consider that f{x,y) represents the evaluation of a physical field 
parametrized by x by a sensor located at a point y. Discarding would model the failure of 
the corresponding sensor. One might still want to use the information provided by the physical 
field parametrized by Xmo which has been selected together with the sensor ymo during the 
greedy procedure. The motivation is that even if the first selected sensor fails, we can still 
use the information from the associated physical field in the construction of the approximate 
separated representation. 

Our idea, inspired from Lemma l2.1l is to choose D = {F'^)f where denotes the pseudo¬ 
inverse. The matrices D and F are both rectangular and in with do = card(A/o). We 

recall that if A € R'^1^'^2^ then A^ € ]^d2xdi is the unique matrix verifying (i) AA^A = A, 
(ii) A^AA^ = Af (hi) (AA^)^ = AAf (iv) (A^A)^ = A^'A. Choosing D = {F'^)f the 
interpolation properties (f3all - (pbl) will not be verified anymore, but since the goal is to find 
an approximation formula and not necessarily an interpolation formula, we can define the 
rectangular EIM approximation by 

(9) 

iGAf mGA/b 

As a numerical illustration, fix E = (1 2 3)^ and consider the function R^ x R 9 (x, y) 1 —>■ 
/(x, y) := cos((E- x)y) € (—1,1). Suppose that d = 8 couples of points {xk,yk), k £ N = 
{l,...,d}, have been selected by the greedy procedure according to (fTa1) - (llbp . Consider any 
pair i,j £ M, i j and set A/q = N \ {i,j}. The relative £^-norm error on 1000 sampling 
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points in ( 0 , 1 )^ x ( 0 , 1 ) is computed by the two following algorithms: (i) using the classical 
EIM approximation with the 6 couples of points {xk,yk), k € A/q, and (ii) using the approx¬ 
imation ([9]) with the 6 points {xk}keAfo ™ x-dimension and the 8 points {yk}k£jV in the 
y-dimension. The relative ^^-norm error has been computed for all A/q = Af\ j}, i,j G AA 
such that i ^ j. In case (i), the maximum, minimum, and mean relative £^-norm errors are 
respectively 6.8 x 10“^, 3.0 x 10“®, and 3.6 x 10“®, whereas in case (ii), these errors are re¬ 
spectively 2.3 X 10“®, 7.6 X 10“^, and 2.4 x 10“®. This example illustrates the fact that using 
the rectangular approximation so as to keep the information from Xi and Xj, i,j G A/" such 
that i ^ j, yields a better accuracy than discarding this information and using the square 
approximation. 

3.2 Application to GEIM |3i| 

Consider a set A” of functions / defined over C M, and a dictionary S of linear forms 
a defined over A. Consider d functions (fi,fd) G A'^ and d linear forms (cri, ,ad) G 
and define the matrix F = (A(/m))i</,me<d' 

Lemma 3.1 (Existence for GEIM) Assume that F is invertible 
the linear combination 

kdif) = ^ Di^rnCriiDfm 

satisfies for all m G {1,... , d}, 

<7m(/) = for all / G A, (11a) 

(y{fm) = cr{Id{fm)), for all a G S. (11b) 

Proof pTa]) implies that amoif) = CTmoildif)) = Y^l<l<d (j2l<m<d DbrnfXmoifrnfj crfif) 
and ([TTb]) implies that a{fmo) = cr{Id{fmo)) = Ei<m<d (Ei<z<d A,mO-z(/mo)) cr{fm), mean¬ 
ing that D = F~'^ can be chosen. □ 

Lemma 3.2 (Stopping criterion for GEIM) Assume that the matrix F^ := {o'i{fm))i<i m<k 

is invertible. If \ak+i{fk+i) - crk+i{Ik{fk+i))\ / 0, then the matrix F'^+^ := {ai{fm))i<i,m<k+i 
is invertible. 

Proof Assume that the matrix A^+^ is not invertible. Then there exists (Ai,..., A^+i) not 
all zero and with A^+i 7 ^ 0, such that Ei<m<fc ^m(Xiifm) = -\k+i(Xi{fk+i), V/ G {1,..., /c + 
1}. We have T.i<i<k Di,moT.i<m<k ^mcriifm) = -Afe+i Ei<«<fc where = 

(^fc)-T definition of and ([TOD, we deduce that Ei<mo<fc Amo/mo = -Afc+i4(/fc+i)- 

Applying fifc+i to both sides of this relation gives ak+i{fk+i) = CFk+i{Ik{fk+i)), yielding the 
desired contradiction. □ 

The generalization of Theorem 11.21 to GEIM, for which any norm on A can be taken, 
follows immediately. 

Let us now return to the setting of sensor failure discussed in Section 13.11 This setting 
nicely fits the GEIM. Suppose that for all u G S, there exists a unique x^ G H such that 
^{f) = fi^a)- In this setting, the functions / represent the parametrized physical field 
depending on x, and cr{f) = f{xa) represents a sensor located at the point x^, that measures 


. Then, setting D = F ^, 
( 10 ) 
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the value of / at this point. Discarding the broken sensor mo, we can produce the following 
rectangular approximation 

E E (12) 

iGAfo m&M 

which should deliver a more accurate representation of / than if the functions fmo were also 
discarded as in the classical (square) EIM context. 
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